rm(list=ls())
library(nomep)
library(coda)

### Set MCMC parameters
BURN<-1000000
MCMC<-5000000
THIN<-1000
VERB<-100000
CHAINS<-8

### Load merged data and drop extraneous variables
load("../data/merged.RData")
covars <- c("big3", "propseats", "ingov", "pro_anti_eu", "emph.europe",
            "country")
cand.raw <- cand.raw[, c(covars, "listnum", "legitcons14", "pirdeucode", 
                         "region")]

### Drop parties with missing data
w.miss <- apply(cand.raw, 1, function (x) any(is.na(x)))
cand.raw <- cand.raw[!w.miss,]
gamma <- gamma[!w.miss,]

### Drop Ireland since it is STV
ireland <- cand.raw$country=="Ireland"
cand.raw <- cand.raw[!ireland,]
gamma <- gamma[!ireland,]

### Build list rank sets and counts for model
cand.raw$np <- NA
cand.raw$Ip <- NA
for (p in unique(cand.raw$pirdeucode)) {
  tmp <- subset(cand.raw, pirdeucode==p)
  for (r in unique(tmp$region)) {
    tmp2 <- subset(tmp, region==r)
    cand.raw$np[cand.raw$pirdeucode==p & cand.raw$region == r] <- sum(tmp2$legitcons14)
    cand.raw$Ip[cand.raw$pirdeucode==p & cand.raw$region == r] <- nrow(tmp2)
  }
}

### Custom fix for list 1642401 NATION.  Sets listnum 2 to 1.
for (i in 1:nrow(cand.raw)) {
  if (cand.raw$pirdeucode[i] == "1642401" & cand.raw$region[i] =="NATION" &
        cand.raw$listnum[i] == 2)
    cand.raw$listnum[i] <- 1
}

### Limit to one entry per list, build covariate matrix, setup Ip,np
cand.raw <- cand.raw[cand.raw$listnum == 1,]
y <- cand.raw[, covars]
y <- cbind(1, y)
colnames(y)[1] <- "(int)"
Ip <- cand.raw$Ip
np <- cand.raw$np

### dummy out countries
country.dummy <- model.matrix(~country, data=y)[,-1]
country.sum <- apply(country.dummy, 2, sum)
colnames(country.dummy) <- sub("country", "", colnames(country.dummy))

### Build final y matrix
y <- y[, c("(int)",
           "emph.europe",
           "big3",
           "ingov",
           "propseats",
           "pro_anti_eu")]
y <- cbind(y, country.dummy)

y.unscaled <- y

### Scale non-dummy variables
y$pro_anti_eu <- scale(y$pro_anti_eu, scale=max(abs(y$pro_anti_eu)))
y$emph.europe <- scale(y$emph.europe, scale=max(abs(y$emph.europe)))

### Fit model
library(parallel)
seed <- 12345
for (i in 1:CHAINS) {
    beta.start <- rnorm(ncol(gamma) * ncol(y), 0, 25)
    mcparallel(plist(as.matrix(y), gamma, Ip, np, burnin=BURN,
                     mcmc=MCMC, thin=THIN, verbose=VERB, c=.09,
                     seed=seed, beta.start=beta.start))
    seed = seed + 1
}
system.time(post <- mccollect())
eval(parse(text=paste("posts <- mcmc.list(", paste("post[[", 1:CHAINS, 
                                                   "]]$beta[,-(1:", ncol(y), ")]",
                                                   sep="", collapse=", "), ")", 
                      sep="")))

save.image("../fits/main.RData")
